This analysis mirrors the earlier one except that is uses the actual areas impacted by PG&E PSPS in 2019, rather than Fire Hazard Severity Zones, with two (2) statewide datasets describing health equity:
Spatial data of PSPS affected areas initially received from PG&E was cross-referenced with reports on their website.
There are some reports on PG&E’s website for which we do not have spatial data about affected areas.
Spatial data from PG&E contained boundaries for 12 separate areas, which - using layer names - were cross-referenced to reports for 7 specific de-energization events (shown on the map below, click for PG&E’s estimates start/end times and number of affected customers).
event1 <- filter(pge_events, event == "June_7-9")
event2 <- filter(pge_events, event == "September_25-27")
event3 <- filter(pge_events, event == "October_5")
event4 <- filter(pge_events, event == "October_9-12")
event5 <- filter(pge_events, event == "October_23-25")
event6 <- filter(pge_events, event == "October_26-November_1")
event7 <- filter(pge_events, event == "November_20-21")
pal2 <- colorQuantile(palette = c("#2172B4","#9DCAE1","#A4DBA1","#008836"), n=4, domain = 0:100)
maphpi <- merge(tracts_sp, {
hpi[,.(ct10, hpi2_pctile)]
})
leaflet() %>%
# addProviderTiles(providers$Esri.WorldStreetMap, group = "Street Map") %>%
# addProviderTiles(providers$Stamen.Terrain, group = "Terrain") %>%
# addProviderTiles(providers$Esri.WorldImagery, group = "Sattelite") %>%
addProviderTiles(providers$Stamen.Toner, group = "B/W") %>%
addPolygons(data = maphpi,
color = "#444444",
weight = 1,
smoothFactor = 0.1,
fillOpacity = 0.5,
fillColor = ~ pal2(maphpi$hpi2_pctile),
stroke = FALSE,
popup = paste0("This location is in the top ",round(maphpi$hpi2_pctile,0),"the percentile of all tracts in the state for the healthy places index"),
group = "HPI"
) %>%
addPolygons(
data = st_zm(event1),
color = "#252525",
weight = 2,
smoothFactor = 0.1,
fillOpacity = 0.7,
fillColor = "#377eb8",
stroke = TRUE,
label = ~event1$event,
popup = paste0("This PSPS event spanned ", round(event1$Duration_hrs,0), " hours \n (beginning ", as.character(event1$start_time), " and ending at ", event1$end_time,"). \n Ultimately PG&E estimates that it affected ",format(event1$customers, big.mark = ',')," customers."),
group = "June 7-9 (blue)"
) %>%
addPolygons(
data = st_zm(event2),
color = "#252525",
weight = 2,
smoothFactor = 0.1,
fillOpacity = 0.7,
fillColor = "#4daf4a",
stroke = TRUE,
label = ~event2$event,
popup = paste0("This PSPS event spanned ", round(event2$Duration_hrs,0), " hours \n (beginning ", as.character(event2$start_time), " and ending at ", event2$end_time,"). \n Ultimately PG&E estimates that it affected ",format(event2$customers, big.mark = ',')," customers."),
group = "Sep 25-27 (green)"
) %>%
addPolygons(
data = st_zm(event3),
color = "#252525",
weight = 2,
smoothFactor = 0.1,
fillOpacity = 0.7,
fillColor = "#984ea3",
stroke = TRUE,
label = ~event3$event,
popup = paste0("This PSPS event spanned ", round(event3$Duration_hrs,0), " hours \n (beginning ", as.character(event3$start_time), " and ending at ", event3$end_time,"). \n Ultimately PG&E estimates that it affected ",format(event3$customers, big.mark = ',')," customers."),
group = "Oct 5 (purple)"
) %>%
addPolygons(
data = st_zm(event4),
color = "#252525",
weight = 2,
smoothFactor = 0.1,
fillOpacity = 0.7,
fillColor = "#ffff33",
stroke = TRUE,
label = ~event4$event,
popup = paste0("This PSPS event spanned ", round(event4$Duration_hrs,0), " hours \n (beginning ", as.character(event4$start_time), " and ending at ", event4$end_time,"). \n Ultimately PG&E estimates that it affected ",format(event4$customers, big.mark = ',')," customers."),
group = "Oct 9-12 (yellow)"
) %>%
addPolygons(
data = st_zm(event5),
color = "#252525",
weight = 2,
smoothFactor = 0.1,
fillOpacity = 0.7,
fillColor = "#a65628",
stroke = TRUE,
label = ~event5$event,
popup = paste0("This PSPS event spanned ", round(event5$Duration_hrs,0), " hours \n (beginning ", as.character(event5$start_time), " and ending at ", event5$end_time,"). \n Ultimately PG&E estimates that it affected ",format(event5$customers, big.mark = ',')," customers."),
group = "Oct 23-25 (brown)"
) %>%
addPolygons(
data = st_zm(event6),
color = "#252525",
weight = 2,
smoothFactor = 0.1,
fillOpacity = 0.7,
fillColor = "#e41a1c",
stroke = TRUE,
label = ~event6$event,
popup = paste0("This PSPS event spanned ", round(event6$Duration_hrs,0), " hours \n (beginning ", as.character(event6$start_time), " and ending at ", event6$end_time,"). \n Ultimately PG&E estimates that it affected ",format(event6$customers, big.mark = ',')," customers."),
group = "Oct 26 - Nov 1 (red)"
) %>%
addPolygons(
data = st_zm(event7),
color = "#252525",
weight = 2,
smoothFactor = 0.1,
fillOpacity = 0.7,
fillColor = "#ff7f00",
stroke = TRUE,
label = ~event7$event,
popup = paste0("This PSPS event spanned ", round(event7$Duration_hrs,0), " hours \n (beginning ", as.character(event7$start_time), " and ending at ", event7$end_time,"). \n Ultimately PG&E estimates that it affected ",format(event7$customers, big.mark = ',')," customers."),
group = "Nov 20-21 (orange)"
) %>%
addLayersControl(
# baseGroups = c("Street Map", "Terrain", "Sattelite", "B/W"),
overlayGroups = c("HPI","June 7-9 (blue)", "Sep 25-27 (green)", "Oct 5 (purple)","Oct 9-12 (yellow)", "Oct 23-25 (brown)" , "Oct 26 - Nov 1 (red)", "Nov 20-21 (orange)"),
options = layersControlOptions(collapsed = FALSE)
)%>% # adds a button to return to the whole state view
addLegend("bottomleft",
pal = pal2,
values = na.omit(maphpi$hpi2_pctile),
labels = c("bottom 25%","","","healthiest 25%"),
opacity = 1,
title = "HPI Percentile") %>%
setView(lat = 39.302394, lng = -122.084003, zoom = 6) %>% #defaults the view of the original map to show the entire state
addEasyButton(easyButton(
icon="fa-globe", title="Zoom to State",
onClick=JS("function(btn, map){ map.setView([37.085206, -119.540085],5); }"))) %>%
hideGroup(c("June 7-9 (blue)", "Sep 25-27 (green)", "Oct 5 (purple)","Oct 9-12 (yellow)", "Oct 23-25 (brown)" , "Oct 26 - Nov 1 (red)", "Nov 20-21 (orange)"))
# adds a button to return to the whole state view
# addLegend("bottomleft",
# pal = pal,
# values = factor(c("Moderate","High","VeryHigh"), levels = c("Moderate","High","VeryHigh")),
# opacity = 1,
# title = "Fire Hazard Severity Zones"
# ##lables = c("X%-X% (most vulnerable)","X%-X%","X%-X%","X%-X%","X%-X% (least"%. In this tract, the "mapTemp$def," is higher than "mapTemp),))##
# )
The census tracts impacted by PSPS can be quickly analyzed using existing Health Equity data from the Healthy Places Index (HPI).
The mean HPI of locations impacted by PSPS tend to be higher (healthier) than those not impacted.
foo <- merge({
master_PSPS_tract[, .(shutoff = ifelse(PSPS_yn == "Yes", 1,0)), by = (GEOID) ] %>% .[, .(PSPS_yn = max(shutoff)), by = (GEOID) ] %>% .[, .(geoid = GEOID, PSPS_yn = ifelse(PSPS_yn ==0,"Not Affected by PSPS", "Affected by PSPS"))]
}, hpi)
foo <- foo[,.(ct10, hpi2score, economic, education, housing, healthcareaccess, neighborhood, pollution, transportation, social, PSPS_yn)] %>%
melt.data.table(id.vars = c("ct10","PSPS_yn"), variable.name = "hpi_component",value.name = "value")
hcboxplot(var = foo$hpi_component, x = foo$value, var2 = foo$PSPS_yn,
outliers = FALSE) %>%
hc_chart(type = "column") %>%
hc_add_theme(hc_theme_smpl()) %>%
hc_xAxis(title = list(text = "HPI Index and Components")) %>%
hc_yAxis(title = list(text = "HPI Score (higher = healthier)")) %>%
hc_title(text = "Healthy Places Index and PSPS-affected areas",
margin = 20, align = "left") %>%
hc_subtitle(text = "data from PG&E and Public Health Alliance of Southern CA",
align = "left", style = list(fontWeight = "bold"))
When we examine the number of times a place has been impacted, it appears that the places impacted the most in 2019 are worse off than those that were impacted once or twice, and in some cases worse than those places not impacted by PSPS.
This result is in part a function of the size of the area and population affected by the different PSPS events (larger events captured more areas that were only impacted by 1 or 2 PSPS events). We would like to examine the idea presented here using duration (hrs) of impact rather than frequency (#), particularly since different events last for different amounts of time. To do this we will need to improve the precision and accuracy of the data from PG&E.
foo <- merge({
master_PSPS_tract[, .(shutoff = ifelse(PSPS_yn == "Yes", 1,0)), by = (GEOID) ] %>%
.[, .(PSPS_yn = max(shutoff), PSPS_count = sum(shutoff)), by = (GEOID) ] %>%
.[, .(geoid = GEOID, PSPS_count = ifelse(PSPS_count >=5, "5 or more",PSPS_count), PSPS_yn = ifelse(PSPS_yn ==1, "Affected by PSPS","Not Affected by PSPS"))]
}, hpi)
foo <- foo[,.(ct10, hpi2score, economic, education, housing, healthcareaccess, neighborhood, pollution, transportation, social, PSPS_yn, PSPS_count)] %>%
melt.data.table(id.vars = c("ct10","PSPS_yn","PSPS_count"), variable.name = "hpi_component",value.name = "value")
hcboxplot(var = foo$hpi_component, x = foo$value, var2 = foo$PSPS_count,
outliers = FALSE) %>%
hc_chart(type = "column") %>%
hc_add_theme(hc_theme_smpl()) %>%
hc_xAxis(title = list(text = "HPI Index and Components")) %>%
hc_yAxis(title = list(text = "HPI Score (higher = healthier)")) %>%
hc_title(text = "Healthy Places Index and PSPS-affected areas",
margin = 20, align = "left") %>%
hc_subtitle(text = "data from PG&E and Public Health Alliance of Southern CA",
align = "left", style = list(fontWeight = "bold"))
We could take a similar look at the data using each of the indicators in the Climate Change and Health Vulnerability Indicators Dataset (note that for CCHVIs high values indicate more vulnerability or worse off conditions).
foo <- merge({
master_PSPS_tract[, .(shutoff = ifelse(PSPS_yn == "Yes", 1,0)), by = (GEOID) ] %>%
.[, .(PSPS_yn = max(shutoff)), by = (GEOID) ] %>%
.[, .(COUNTYFI_1 = GEOID, PSPS_yn = ifelse(PSPS_yn ==0,"Not Affected by PSPS", "Affected by PSPS"))]
}, cchvi_tracts[race == "Total" & latest =="Y"], by = "COUNTYFI_1")
foo <- foo[,.(ct10, est, ind_strt, PSPS_yn)] #%>%
#melt.data.table(id.vars = c("ct10","PSPS_yn"), variable.name = "hpi_component",value.name = "value")
hcboxplot(var = foo$ind_strt, x = foo$est, var2 = foo$PSPS_yn,
outliers = FALSE) %>%
hc_chart(type = "column") %>%
hc_add_theme(hc_theme_smpl()) %>%
hc_xAxis(title = list(text = "CCHIV Indicators")) %>%
hc_yAxis(title = list(text = "CCHVI Values (higher = more vulnerable)")) %>%
hc_title(text = "Climate Change & Health Vulnerability Indicators and PSPS-affected areas",
margin = 20, align = "left") %>%
hc_subtitle(text = "data from PG&E and CDPH",
align = "left", style = list(fontWeight = "bold"))
Each of the plots above are boxplots, which describe the distribution of each indicator in the different areas. The center line in those plots represent the median for the group, and the whiskers are extremes. Moving forward, specific indicators (by event) could be investigated in a more rigorous way using data from the American Community Survey (ACS) directly to compare the estimated indicator values for different areas and to provide confidence intervals, which allow for statistical testing of differences between PSPS-affected areas and non-affected areas.
Below is an example for the % of people over 65 in the different areas by event. Similar comparisons have also been done for the number of adults living with a disability (and below the poverty line).
old <- tidycensus::get_acs(geography = "tract",
variables = c(numerator = "B01001_020",
numerator = "B01001_021",
numerator = "B01001_022",
numerator = "B01001_023",
numerator = "B01001_024",
numerator = "B01001_025",
numerator = "B01001_044",
numerator = "B01001_045",
numerator = "B01001_046",
numerator = "B01001_047",
numerator = "B01001_048",
numerator = "B01001_049",
denominator = "B01001_001"),
state = "CA") %>% data.table()
old <- old[,.(estimate = sum(estimate),
moe = tidycensus::moe_sum(estimate = estimate, moe = moe)),
by= .(GEOID, NAME, variable)]
try1 <- merge(old, {
readRDS("//mnt/projects/ohe/PSPS/PSPS_yn_master_tract.RDS")
}, allow.cartesian = T)
# try1 <- na.omit(try1)
try1 <- try1[, se:= moe/1.645]
sums <- try1[,.(sum = sum(estimate, na.rm=T),
se_sum = (sum(se^2, na.rm=T)^0.5)), by = .(variable, PSPS_yn, event)]
totals <- sums %>% dcast.data.table(formula = PSPS_yn + event ~ variable, value.var = "sum")
sum_ses <- sums %>% dcast.data.table(formula = PSPS_yn + event ~ variable, value.var = "se_sum") %>% .[,.(PSPS_yn, event, SE_num = numerator, SE_pop = denominator)]
final <- merge(totals, sum_ses)
final <- final[, `:=` (proportion = numerator/denominator,
se_prop = (1/denominator)*(SE_num^2 - (((numerator/denominator)^2)*SE_pop^2))^0.5)]
final <- final[, `:=` (LL = proportion - 1.96*se_prop,
UL = proportion + 1.96*se_prop)]
final <- final[, event := factor(event, levels = c("June_7-9", "September_25-27","October_5","October_9-12","October_23-25","October_26-November_1","November_20-21"))] %>% setorder(event)
final <- final[, PSPS_yn := ifelse(PSPS_yn == "No","Not Affected by PSPS", "Affected by PSPS")]
hc1 <- highchart() %>%
hc_chart(type = "bar") %>%
hc_add_series(final, "point", hcaes(y = 100*proportion, x = event, group = PSPS_yn))%>%
hc_add_series(final, "errorbar", hcaes(low = 100*LL, high = 100*UL, x = event, group = PSPS_yn)) %>%
hc_add_theme(hc_theme_smpl()) %>%
hc_xAxis(title = list(text = "PSPS Event"),
type = 'category' ,
categories = c("June 7-9 (14hrs, 22k)", "Sep 25-27 (26hrs, 75k)","Oct 5 (18hrs, 11k) ","Oct 9-12 (90hrs, 729k)","Oct 23-25 (52hrs, 177k)","Oct 26-Nov 1 (129hrs, 941k)","Nov 20-21 (42hrs, 49k)")
) %>%
hc_yAxis(title = list(text = "% over 65 (with 95% CI)")) %>%
hc_title(text = "% of population over 65",
margin = 20, align = "left") %>%
hc_subtitle(text = "Tracts Impacted by PSPS vs Not Impacted",
align = "left", style = list(fontWeight = "bold"))
hc1
This Table summarizes much of the data we currently have from PG&E on their events as well a few of the comparisons on indicators of vulnerability for areas impacted and not impacted in each event.
| Jun 7-9 | Sep 25-27 | Oct 5 | Oct 9-12 | Oct 23-25 | Oct 26-Nov 1 | Nov 20-21 | ||
|---|---|---|---|---|---|---|---|---|
| start time | 6/8/2019 6:00 | 9/23/2019 17:06 | 10/5/2019 22:00 | 10/9/2019 0:09 | 10/23/2019 13:54 | 10/26/2019 16:00 | 11/20/2019 3:38 | |
| end time | 6/8/2019 20:00 | 9/24/2019 18:40 | 10/6/2019 16:00 | 10/12/2019 17:41 | 10/25/2019 18:20 | 11/1/2019 1:20 | 11/21/2019 21:55 | |
| Duration_hrs | 14 | 25.57 | 18 | 89.53 | 52.43 | 129.33 | 42.28 | |
| customers | 22,327 | 75,385 | 11,304 | 728,980 | 176,620 | 941,217 | 49,085 | |
| residential | 19,500 | 35,501 | 9,981 | 636,355 | 832,314 | |||
| comm_indust | 2,565 | 3,856 | 1,250 | 81,318 | 98,514 | |||
| med_baseline | 1,589 | 4,451 | 718 | 30,026 | 7,823 | 34,618 | 2,456 | |
| other | 262 | 616 | 73 | 11,307 | 10,389 | |||
| % Over 65 | ||||||||
| Affected by PSPS | 22.8% (22.1, 23.5) | 23.6% (22.8, 24.4) | 24.3% (23, 25.7) | 17.9% (17.7, 18) | 24.5% (24, 25) | 18.7% (18.6, 18.9) | 24.2% (23.3, 25.1) | |
| Not Affected | 13.2% (13.1, 13.2) | 13.2% (13.1, 13.2) | 13.2% (13.2, 13.2) | 12.9% (12.8, 12.9) | 13.1% (13, 13.1) | 12.7% (12.7, 12.8) | 13.2% (13.1, 13.2) | |
| % Disabled | ||||||||
| Affected by PSPS | 11.1% (10.5, 11.8) | 19.6% (18.7, 20.5) | 23.2% (21.5, 24.9) | 11.9% (11.8, 12.1) | 15.9% (15.4, 16.4) | 11.6% (11.5, 11.8) | 16.9% (15.8, 17.9) | |
| Not Affected | 9.9% (9.9, 10) | 9.9% (9.8, 9.9) | 9.9% (9.9, 9.9) | 9.8% (9.7, 9.8) | 9.8% (9.8, 9.9) | 9.8% (9.7, 9.8) | 9.9% (9.9, 9.9) | |
| % Disabled and Below Poverty Line | ||||||||
| Affected by PSPS | 1.5% (1.3, 1.8) | 3.8% (3.4, 4.2) | 4.4% (3.6, 5.2) | 2.1% (2, 2.2) | 2.5% (2.3, 2.8) | 1.9% (1.8, 2) | 3.2% (2.7, 3.7) | |
| Not Affected | 2% (2, 2) | 2% (1.9, 2) | 2% (2, 2) | 2% (1.9, 2) | 2% (1.9, 2) | 2% (2, 2) | 2% (2, 2) |